program VMD_LIN

  use fft_wrap
  use mpi_wrap
  use efield
  use collisions

!  use evolve

  use dimensions

  use parameters, only : restart, nt, dt, iDiag, iData, vPhase, runname

  implicit none


  !-- major arrays ---------------------

  complex, allocatable  :: dfu(:, :, :, :, :)   ! distribution function
  complex, allocatable  :: uxMat(:, :, :, :)    ! velocity, x-component
  complex, allocatable  :: uyMat(:, :, :, :)    ! velocity, y-component, unused
  complex, allocatable  :: Ex2D(:, :, :)        ! electric field, x-component 
  complex, allocatable  :: Ey2D(:, :, :)        ! electric field, y-component 
  complex, allocatable  :: Ez2D(:, :, :)        ! electric field, z-component

  real, allocatable     :: filter(:,:,:)        ! conical filter

  real, allocatable     :: Fbgk(:,:)            ! f_bgk
  complex, allocatable  :: ftDVbgk(:,:)         ! F.T. in v of (df_bgk / dv)
  real, allocatable     :: Ebgk(:)              ! E_bgk
  real, allocatable     :: Ubgk(:)              ! longitudinal velocity



  !-- variables  ----------------------- 

  real       :: t                               ! time
  integer    :: it,iz                           ! iteration count
  integer    :: iargc                           ! argument count


!-------------------------------------------------------------------------------------

  allocate ( dfu(jv, jz, jy, jx, jDir) )
  allocate ( uxMat(jz, jy, jx, jDir) )
  allocate ( uyMat(jz, jy, jx, jDir) )
  allocate ( Ex2D(jz, jy, jx) )
  allocate ( Ey2D(jz, jy, jx) )
  allocate ( Ez2D(jz, jy, jx) )

  allocate ( filter(jz,jy,jx) )

  allocate ( Fbgk(jv,jz) )
  allocate ( ftDVbgk(jv,jz) )
  allocate ( Ebgk(jz) )
  allocate ( Ubgk(jz) )


      !--  reading the run_name from the command line

      if(iargc().eq.0) then
         call getarg(0,runname)
         write(*,*)      'Usage: ',trim(runname),' <run name>'
         call mpi_abort(MPI_COMM_WORLD, -1, mpi_err)
         stop
      end if
      call getarg(1,runname)


      !-- set parameters and coordinate space --

      call mpi_wrap_init()

      call read_input()

      call fft_init(dfu)

      call set_space_vectors()

      call efield_init()

      !-- restore from restart file or create initial conditions --

      call maxwell(Fbgk)

      call io_read_bgk(Fbgk, dfu)

      !call bgk_test_equil(Fbgk, DVbgk, Ebgk)
      !call mpi_wrap_exit()
      !return

      call efield_bgk(Fbgk, Ebgk)
      call derivV_bgk(Fbgk, ftDVbgk)
      call sumUz_bgk(Fbgk,  Ubgk)

      !-- print out Uz0(z) for debugging --
      !open(93, file= 'z_rho_Uz_E.txt')
      !do iz = 1,jz
      !write (93,*) zvec(iz), sum(Fbgk(:,iz))*dv,  Ubgk(iz), Ebgk(iz)
      !end do
      !close(93)
      !return

      Fbgk  = Fbgk / (jDir*npDir)
      ftDVbgk = ftDVbgk /(jDir*npDir)

      call set_ic_lin(dfu, uxMat, Fbgk, it)

      t  = it*dt

      call set_filter(filter)

      call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)
      call diagnostics_lin( dfu, Ex2D, Ez2D, uxMat, t)

      call write_bin(dfu, uxMat, uyMat, it, it/iData)

      write(msg,*) char(10), ' Begin temporal evolution', char(10)
      call mpi_wrap_msg


!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------

  call evolve_F_collisions(dfu, dt/2)

  do while (it .le. nt) 

    call dealias_5D(dfu)
    call dealias_4D(uxMat)
    call dealias_4D(uyMat)


    call evolve_lin_F_long_advect( dfu, dt/2)
    call evolve_lin_F_trans_advect( dfu, Fbgk, uxMat, dt/2)

    !call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)
    !call evolve_lin_F_accel( dfu, Ebgk, ftDVbgk, Ez2D, dt)
    !call evolve_lin_U (uxMat, Ex2D, Ubgk, dt)

    !-- three commented out lines above are replaced by:
    !
    call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)
    call evolve_lin_U (uxMat, Ex2D, Ubgk, dt/2)
    call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)
    call evolve_lin_F_accel( dfu, Ebgk, ftDVbgk, Ez2D, dt)
    call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)
    call evolve_lin_U (uxMat, Ex2D, Ubgk, dt/2)

 
    call filter_U_lin( uxMat, filter )

    call evolve_lin_F_trans_advect( dfu, Fbgk, uxMat, dt/2)
    call evolve_lin_F_long_advect( dfu, dt/2)

    t  = t + dt
    it = it + 1
   
    !-- diagnostics --

    if ( (it/iDiag)*iDiag.eq.it ) then      

      call evolve_F_collisions(dfu, dt/2)

      call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)

      call diagnostics_lin( dfu, Ex2D, Ez2D, uxMat, t)
      !write(msg, "(' Diagnostcs, t =' e12.6)") t
      !call mpi_wrap_msg
  
      !-- output binary data for restart and post-processing --

      if ( (it/iData)*iData.eq.it ) then

         call write_bin(dfu, uxMat, uyMat, it, it/iData)

         write(msg, "(' Output ' i4.4)") it/iData
         call mpi_wrap_msg
       
      end if

      call evolve_F_collisions(dfu, dt/2)

    else

      call evolve_F_collisions(dfu, dt)

    end if

  end do      ! END TEMPORAL LOOP
 

  write(msg,*) char(10), ' End temporal evolution', char(10)
  call mpi_wrap_msg

  !-- write some specrtral spectrum quantities --
       
  call mpi_wrap_exit()
    
  stop;

end


!-------------------------------------------------------------------------------------
subroutine maxwell(Fbgk)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname, Vmax, pi, vPhase, iWaveFrame

  implicit none

  real, intent(out)       :: Fbgk(jv,jz)                     ! f_bgk
  !real, intent(out)       :: DVbgk(jv,jz)                    ! df_bgk / dv
  !real, intent(out)       :: Ebgk(jz)                        ! E_bgk

  real                    :: f0
  integer                 :: iv

  !Ebgk(:)=0.0

  ! reminder:  vvec(iv) = -vPhase*iWaveFrame - Vmax + dv*(iv-1)
  ! reminder:  f0vec(:) = ( 1/sqrt(2*pi) )*exp(-0.5*(vvec(:) + vPhase*iWaveFrame)**2)

  do iv = 1,jv
     f0 = ( 1/sqrt(2*pi) )*exp(-0.5*(dv*(iv-1) - Vmax)**2)
     Fbgk(iv,:) = f0
     !DVbgk(iv,:) =  -(dv*(iv-1) - Vmax)*f0
  end do

end subroutine maxwell

!-------------------------------------------------------------------------------------
